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ABSTRACT 

Pulsars are famous for their rotational stability. Most of them steadily spin down and dis¬ 
play a highly repetitive pulse shape. But some pulsars experience timing irregularities such 
as nulling, intermittency, mode changing and timing noise. As changes in the pulse shape are 
often correlated with timing irregularities, precession is a possible cause of these phenomena. 
Whereas pulsar magnetospheres are filled with plasma, most pulsar precession studies were 
carried out within the vacuum approximation and neglected the effects of magnetospheric 
currents and charges. Recent numerical simulations of plasma-filled pulsar magnetospheres 
provide us with a detailed quantitative description of magnetospheric torques exerted on the 
pulsar surface. In this paper, we present the study of neutron star evolution using these new 
torque expressions. We show that they lead to (1) much slower long-term evolution of pul¬ 
sar parameters and (2) much less extreme solutions for these parameters than the vacuum 
magnetosphere models. To facilitate the interpretation of observed pulsar timing residuals, we 
derive an analytic model that (1) describes the time evolution of non-spherical pulsars and (2) 
translates the observed pulsar timing residuals into the geometrical parameters of the pulsar. 
We apply this model to two pulsars with very different temporal behaviours. For the pulsar 
B1828-11, we demonstrate that the timing residual curves allow two pulsar geometries: one 
with stellar deformation pointing along the magnetic axis and one along the rotational axis. 
For the Crab pulsar, we use the model show that the recent observation of its magnetic and 
rotational axes moving away from each other can be explained by precession. 
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1 INTRODUCTION 

Studies of pulsar evolution are often performed under the 
assumption of a spherical neutron star in vacuum. In ad¬ 
dition to such idealized studies (Deutsch 1955), it is now 
possible to perform more realistic numerical simulations in 
the limit of a plasma-filled magnetosphere (Contopoulos et al. 
1999; Spitkovsky 2006; Kalapotharakos & Contopoulos 2009; 
Tchekhovskoy et al. 2013), and more recently in the full kinetic 
limit (Philippov & Spitkovsky 2014; Chen & Beloborodov 2014; 
Cerutti et al. 2015; Philippov et al. 2015; Belyaev 2015). These 
models predict the evolution of the pulsar rotation period P and the 
inclination angle a that the magnetic axis makes with the rotational 
axes. For a spherically-symmetric star, the pulsar period slowly in¬ 
creases and the inclination angle gradually decreases towards 0° 
(Philippov et al. 2014; but see Beskin et al. 1993). Flowever, devi¬ 
ations in the pulsar pulse arrival times from the expectations of the 
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above smooth time-evolution, e.g., pulsar timing noise and inter¬ 
mittency, indicate that pulsar spindown may experience quite dif¬ 
ferent, often periodical changes. 

One way to introduce an oscillatory behaviour into the spin- 
down is to consider a non-spherical star. In the absence of a phys¬ 
ically motivated model of a plasma-filled magnetosphere, previous 
studies (Goldreich 1970; Melatos 2000) considered pulsar evolu¬ 
tion in the vacuum approximation. In this paper, we investigate 
the evolution of non-spherical pulsars whose magnetospheres are 
filled with plasma. For this, we use recent results of magnetohy¬ 
drodynamic (MHD) simulations of plasma-filled magnetospheres 
(Philippov et al. 2014). 

The paper is organized as follows. We start with the discussion 
of the vacuum approximation for pulsar spindown in Section 2. We 
discuss the parametrisation of magnetospheric torques in Section 3. 
The solutions for spherical pulsars are presented in Section 4. The 
effects of stellar non-sphericity are discussed in Section 5. We dis¬ 
cuss the implications of our results in Section 6 and conclude in 
Section 7. 
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2 VACUUM SOLUTION 


A pulsar rotating in vacuum emits magnetodipole radiation that car¬ 
ries away the angular momentum of the neutron star. The resulting 
torque acting on the star could he found hy integrating the stresses 
exerted at its surface, 


K,= 


f 




( 1 ) 


where the summation over repeated indices is implied, T*/ is the 
stress tensor (in our case, we have only electromagnetic forces act¬ 
ing on the star, so it is Maxwell’s stress tensor), is the Levi- 
Civita symbol, r is the position of surface element with an area 
element d5, and n; is the normal vector. 

If we know electric and magnetic field components near the 
stellar surface, it is straightforward to compute the torque K from 
eq. (1). Deutsch (1955) computed vacuum electromagnetic fields 
near a rotating spherical star endowed with a dipolar magnetic 
field. Using these expressions, Michel & Goldwire (1970) derived 
the torques acting on the pulsar in vacuum. If pulsar’s angular fre¬ 
quency vector n is oriented along the z-axis and the pulsar mag¬ 
netic moment vector// lies in the x-z plane, the torque components 
in the vacuum limit are: 


where 
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2 
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( 2 ) 

(3) 

(4) 


and a is the inclination angle, i.e., the angle between il and//. 

The torque component due to eq. (2) slows down stellar ro¬ 
tation. The torque component given by eq. (3) points towards //. 
Under the action of this torque component, the inclination angle 
a evolves toward the alignment, and the star approaches the state 
of minimum angular momentum and energy loss (Philippov et al. 
2014). 

Though the vacuum magnetosphere model predicts the spin- 
down rate close to the observed value, it has several important 
shortcomings. The magnetospheres of real pulsars are filled with 
plasma that inevitably modifies the structure of the magnetic fields 
via the charges and currents that it carries (Goldreich & Julian 
1969). We discuss these plasma effects below. 


3 MAGNETOSPHERIC TORQUES 

Unfortunately, there is no simple analytical solution for mag- 
netospheric torques in the plasma-filled magnetosphere case. 
Thus, we will rely on the recent advances in numerical sim¬ 
ulations of oblique pulsar magnetospheres (Spitkovsky 2006; 
Kalapotharakos & Contopoulos 2009; Tchekhovskoy et al. 2013; 
Philippov et al. 2015; Tchekhovskoy et al. 2015) that allow us to 
quantitatively study the influence of magnetospheric charges and 
currents on the structure of the magnetosphere as well as on the 
pulsar spindown and alignment. 

Philippov et al. (2014) analyzed the results of force-free and 
MHD simulations and came up with the following parametriza- 
tion of magnetospheric torques for plasma-filled magnetospheres, 
to which for simplicity we will refer to below as simply MHD mag- 


Table 1. Spindown parameters of different magnetospheric models, ko, k[ 
and k 2 are defined in equations (5) - (6), and G is referred to the anomalous 
torque and is defined in equation (7). 


Model 

kg 

ki 

h 

ki 

Vacuum 

0 

2/3 

2/3 

0.3 

MHD/force-free 

1 

1 

1 

~ 0.1 


netospheres: 

A “ ~^aligReddo Sttl O'), (5) 

Aj. = A ;2 Aaiigned sin a cos a, (6) 


where kg, ki and ^2 are factors of order unity. They weakly depend 
on RflRix, the ratio of the stellar radius R, to the light cylinder ra¬ 
dius Ac = cjO.. Typical parameters for vacuum and MHD models 
are listed in Table 1 (Philippov et al. 2014). 

In addition to torque components (5) and (6), we need to know 
the anomalous torque, or the torque component directed along 12 x 
p, which could be written as 


A, 


^ hK-idgRRi-— sin a cos a, 

\IK 


(7) 


where kg = const. Melatos (2000) estimated kg = 0.3. In numeri¬ 
cal simulations, anomalous torque is hard to measure as it diverges 
with J? —> 0. We estimate it in our MHD model to be kg ~ 0.1 
(please see Table 1). As we will show below, the timing residu¬ 
als do not depend on the strength of the anomalous torque and the 
value of kg (see Section 5.3). 

The most important difference between vacuum and MHD 
models is in the value of parameter kg, which describes the energy 
losses of an aligned rotator. In vacuum, an aligned rotator does not 
spin down, i.e. or = 0, 12 0 is a solution to evolution equations. 

This is impossible in the MHD case, because the spindown lumi¬ 
nosity of an aligned rotator is non-zero. 


4 EVOLUTION OF SPHERICAL STARS 


The torque given by equations (5) - (6) forces the star to change its 
angular velocity as 


dO 

dt 


( 8 ) 


where / is the stellar moment of inertia. The torque component A 
acts in the direction opposite to 12 and reduces the magnitude of 
angular velocity without changing its direction. On the other hand, 
A and Ky act perpendicular to 12 and change its direction with¬ 
out changing its magnitude. A moves angular velocity towards 
the magnetic moment p, so it changes the inclination angle a. The 
anomalous torque component, A, causes the star to precess on the 
timescale of 


InlRd 

p-Clkg 


RcPkg ’ 


(9) 


where the prefactor 


^ = kg + ki sin^ a 


( 10 ) 


reflects the inclination-dependence of pulsar radiation losses (see 
eq. 5). As we will see below, the anomalous torque does not affect 
the long-term evolution of spherical stars and has only a limited 
effect on the evolution of non-spherical stars. 
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Figure 1. Evolution of a spherical pulsar with initial angle cro = 60 degrees and angular frequency Qo. Time on the upper horizontal axis is measured in units 
of characteristic pulsar spin down time t = [Left panel:] The evolution of angular velocity of the neutron star. Vacuum pulsars evolve toward a 

constant rotational frequency, £lo cos oq, while MHD pulsars follow power-law dependence, H oc for t:^T. [Right panel:] The evolution of the obliquity 
angle. In the vacuum case, the evolution is exponentially fast and leads to a very sharp decrease in the inclination angle. In the MHD case, the evolution is a 
power-law in time, a oc 


A spherical star evolves according to, 
dfi 
df 
„ dor 

IQ .— = -K^. 
df 


( 11 ) 

( 12 ) 


The solution of pulsar evolution equations could be simplified, be¬ 
cause the system of equations (5)-(6) and (11)-(12) has an integral 


n 






- 




sin*^“ oq 


ap j 


i/h 


(13) 


\ sin*" a 

which in the case of vacuum losses (kg = 0, ki = k 2 - 2/3, see 
Tab. 1) reduces to 


n cos a = no cos ffo. 


(14) 


and in the MHD case (feo = A;i = ^2 = 1) to 
cos^ a cos^ ag 

i2- = i2o-. 

sin or sin 0*0 


(15) 


Most of the works on the pulsar evolution assume magnetodipole 
(vacuum) radiation and ignore the evolution of the obliquity an¬ 
gle, a (Popov & Turolla 2012; Igoshev & Popov 2013). However, 
as one can show from equations (11)-(12), the inclination angle 
evolves on the same timescale. 


IQ. _ Ic’ _ 

^aligned fl^Q^ P 


(16) 


as the angular velocity, where f is given by eq. (10). Here P and P 
are the pulsar period and period time derivative, respectively. 

As pulsars spin down, the obliquity angle evolves and this af¬ 
fects the spindown. Therefore, it is crucial to account for the change 
in a. For instance, if we do not account for the alignment effect in 
the vacuum case, i.e. postulate a = constant, a pulsar starting out 
with n = Do and a = ag evolves towards Henai = 0. However, the 
correct solution of eqs. (11)-(12) with torques given by eqs. (5)-(6) 
yields Hfinai = flocosoro, which is comparable to flo for represen¬ 
tative ag ~ 60°. 

Though both vacuum and MHD models evolve toward the 
aligned configuration, they approach it through quite different 


paths. The inclination angle of a vacuum pulsar evolves according 
to 


sin a = sin ag exp (-f/, (17) 

where = 1.5 tcos“^ ao, and r is pulsar spindown timescale 
given by eq. (16). Thus, the vacuum pulsar evolves to the aligned 
configuration exponentially fast, without a significant slowdown of 
its rotation. For an MHD pulsar, we have the following implicit 
solution. 


-5— -I- log(sin a) - 

2 sin a 


t 


' align 



-I- log(sin Q-o), 


(18) 


where = t sin^ agj cos"^ ag, which asymptotes at late times to 
a power-law, a oc (Philippov et al. 2014). 

Figure 1 illustrates the evolution under vacuum and MHD 
torques of a spherical pulsar with an initial angle of ag = 60 de¬ 
grees and Crab-like parameters, P - 0.033 s, 6 = 3.78 x 10'^ G, 
T = 7.4 X 10^ years. One can see that for a vacuum pulsar the incli¬ 
nation angle becomes exponentially small just after t > which 
equals to 6 t for Crab’s parameters. On the other hand, the obliq¬ 
uity angle in the MHD case changes only by a factor of few for 
the same parameters. Note also that MHD pulsars lose angular mo¬ 
mentum more efficiently than vacuum ones. Instead of pulsar rota¬ 
tion asymptoting to Oflnai = Qg cos ag as for vacuum pulsars, MHD 
pulsars lose a substantial amount of angular momentum during its 
lifetime, and Q continuously decreases, asymptotically approach¬ 
ing zero: D oc Characteristic timescale of pulsar spindown is 
around lO'^-lO^ years (i.e. the timescale over which D changes by 
a factor of 2). 


5 EVOLUTION OF NON-SPHERICAL STARS 

The rotation of pulsars steadily slows down, and the radio pulse 
profile is mostly stable in time. The deviation from smooth spin- 
down is often less than one part in 10", but at this accuracy pulsars 
experience rather random irregularities, which are referred to as 
timing noise. 
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Figure 3. The time evolution of a non-spherical pulsar with the initial inclination angle ag = 60 deg, angles Bq = 60 deg and = 1 deg (see Fig. 2 for the 
angle definitions), ellipticity e = 1.66 X 10“'^, and angular frequency CIq. [Left panel:] The angular velocity O. vs time. Despite the star is non-spherical, D 
evolves in a very similar way to the spherical case. [Right panel:] The inclination angle a vs time. The stellar non-sphericity causes a to oscillate at late time, 
t > lOr in both vacuum and MHD cases. Despite that, in a time-average sense the inclination angle undergoes a secular decrease. Eventually, the inclination 
angle stabilises. The stabilisation in MHD pulsars occurs at a much later time than in vacuum pulsars. 


63 


Figure 2. In this work, we use the following coordinate system: ei , 82 , 
82 denote the principal axes of the star, which is characterised by a stellar 
magnetic moment // and angular frequency il. Pulsar geometry is defined 
by several angles: x between // and 63 , 6 between 42 and 82 , and a (not 
shown) between // and 42. During free precession, 42 rotates around ^3 and 
the angle between the 8 [ and the projection of 42 on 8[-82 plane increases 
at a constant rate tjp. 


Lyne et al. (2010) presented timing residuals (deviations from 
steady spindown) of 17 pulsars that showed a quasi-periodic be¬ 
haviour. The characteristic period of these variations is a few 
hundred days, hence prolonged observations are needed to inves¬ 
tigate this behaviour. Interestingly, a lot of pulsars with timing 
noise exhibit variations in the pulse profile that are correlated with 
variations in the period. This makes a change in geometry (e.g., 
precession-induced changes in the inclination angle) a possible ex¬ 


planation of this phenomenon. Below we focus on the two cases 
of pulsar B1828-11 and the Crab pulsar. In this section, we set 
R = 10^ cm and consider the case of a plasma-filled magneto¬ 
sphere, i.e., set kg = ki = k 2 = 1 . 

Previously, we discussed a spherical star with an isotropic ten¬ 
sor of inertia. Generally, a rigid star has a biaxial inertia tensor that 
is diagonalised in the basis of its principal axes. Let us denote the 
pulsar principal axes as ej, 82 , 82 and corresponding moments of 
inertia as Ii = I, p = /(I + S 12 ) and p = /(I -t- £ 13 ), respectively. 
Here, £12 and £13 characterise the degree of stellar non-sphericity. 

We will start with the case of a biaxial star, i.e., we will set 
£12 = 0, ei 2 = e t 0. Without the loss of generality, we choose e, 
to lie in the 82 -fJ plane. As seen in Fig. 2, the problem geometry is 
now set by three angles: angle x between /t and 82 , angle 9 between 
41 and £3 and angle a between 41 and //. Note that the values of 
these angles are not completely arbitrary but are constrained by 
|0 < a < |6 + x\- Here and below, indices (1, 2, 3) refer to 

the coordinate system set by the principal axes and (x, y, z) to the 
coordinate system described in Section 2. 

For a non-spherical pulsar, instead of equations (11)-(12) we 
need to solve Euler’s equations of motion of a rigid body: 

L, = eij.LjQ., + K„ (19) 

where summation over repeated indices is implied, Kj and 41, are 
the components of K and 41 in the (ei, 82 , £ 3 ) coordinate system 
and Lj = /,42, . 

To compute Kj, we need to make a coordinate transformation 
from (x, y, z) to (ei, 82 , £ 3 ). This transformation could be found in 
Melatos (2000), and we do not give it here. Rewriting equation (19) 
in the new coordinates, we obtain: 


-t- £^22^3 “ K[ll, 

(20) 

^^2 ~ £^ 23 ^ 2 ] = 7 ^ 2 /^! 

(21) 

(1 + £)i 23 = j /. 

(22) 


We solve these equations numerically using the three angles 
00 , To and bo and the initial period Pq as the initial conditions. We 
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adopt the following fiducial parameters for a neutron star: mass 
M - IAMq, radius = 10 km, dipole moment yu = Q.5BS?, and 
B = 10*2 Q 

To see how a small stellar non-sphericity can affect the evo¬ 
lution of a neutron star, let us again consider a Crab-like pulsar, 
tto = 60°, P = 0.033 s, Bi 2 = 3.78, t = 7.4x 10^ years. Fig. 3 shows 
how this pulsar evolves for 9o = 60°,xo = 1° and e = 1.66 x lO^*^. 
The right panel in Fig. 3 shows that even such a seemingly small 
ellipticity causes the inclination angle a to vary with a substantial 
amplitude. These variations are due to the free precession of the 
star, at the period set by the degree of stellar non-sphericity: 

Tprec = 2;TTprec = p (23) 

which for our choice of parameters is Tprec = 6300 years. As we 
will see below, the standard expression (23) applies only in the 
limit of angular frequency pointing along the symmetry axis of the 
neutron star, ej, i.e., for cos 6o ~ 1 [see equation (48), which gen¬ 
eralises eq. (23) to other values of 0o]. 

In vacuum, the rotation of a spherical star asymptotes to a non¬ 
zero value: this is because pulsars quickly evolve toward the align¬ 
ment a = 0 and their spindown torque vanishes. In contrast, the 
left panel in Fig. 3 shows that non-spherical vacuum pulsar rota¬ 
tion slows down all the way to zero. This is because non-sphericity 
effects cause the inclination angle to stabilise at a finite value, e.g., 
at = 0.35 deg for our choice of pulsar parameters, which pre¬ 
vents the spindown torque from vanishing (see eq. 2). Similarly, for 
MHD pulsars, stellar non-sphericity leads to the inclination angle 
stabilisation at a finite value, e.g., at = 3.46° for our choice of 
pulsar parameters, as seen in the right panel of Fig. 3. The value of 
the stabilisation angle depends on the strength of anomalous torque 
and the k-^ factor (see Tab. 1) and can substantially dilfer between 
vacuum and MHD models; a study of this dependence is beyond 
the scope of this work. 

After stabilisation, the evolution of the angular velocity of the 
neutron star can be approximated by 

/Q\3/;2 

(£2) =- 

where brackets (...) denote the averaging over a precession period. 
Equation (24) has a solution 

1 I _ t 

(£2)2 £22^^ Tevol 

where Teyoi is the characteristic timescale of evolution and £2s,ai, is 
the value of angular frequency at stabilisation. The timescale Teyoi 
is determined by the stabilisation angle and is different in vacuum 
and MHD cases: 

2'era? = 1 -St sin“2 .r, (26) 

= ■'■(1 + astab)“‘ ~ T. (27) 

Equation (25) implies a plateau in (£2) at early times post¬ 
stabilisation (see the left panel of Fig. 3) and a power-law evolution 
at late times, t > Teyoi. For the case shown in Figure 3 the analytical 
estimate (26) gives t^^^ = 5.21 x 10* years, in agreement with the 
end of plateau. 

Figure 3 shows that until the inclination angle stabilises, the 
solution for a non-spherical pulsar is similar to the solution for 
a spherically-symmetric neutron star, except for small-amplitude 
short timescale oscillations super-imposed on top of the secular 
alignment trend. This allows us make use of eqs. (17) and (18) 
to obtain a simple estimate of the time at which the stabilisation 


sets m: 


WAC 

^stab 


= 1.5t 


loglsinffo/sinttstab) 


few X T, 


AlHD ‘ 


(28) 

(29) 


Here we assumed the stabilisation angle to be small and neglected 
the logarithmic factor in eq. (28). Importantly, the stabilisation time 
for MHD pulsars is much larger than for vacuum ones: 


.MHD 

^stab 

fVAC 

^stab 


1 

few X 0-2 


» 1 , 


(30) 


but even for the vacuum pulsar this is a large timescale: for the 
initial conditions used in Fig. 3, ~ 10^ years. For MHD pulsars 

it is > 10* years (for ctstab :£ 0.03). We thus expect that pulsars cross 
the death line before their inclination angle has a chance to reach a 
small value and become stabilised. 


5.1 Neutron Star Ellipticities 

Due to the uncertainties in the neutron star equation of state 
(Glendenning 1997), it is hard to estimate how much the pulsar 
crust is deformed and therefore to theoretically compute the value 
of ellipticity. Horowitz & Kadau (2009) performed multi-million 
ion molecular dynamics simulations showing that the neutron star 
crust can support ellipticities up to = 4 x 10“*. This ellipticity 
corresponds to a precession time of a few days for a one-second 
pulsar (see eq. 23). 

There are two main contributors to the neutron star ellipticity: 
fast rotation and strong magnetic fields in the interior. The non¬ 
sphericity caused by rotation can be estimated as the ratio between 
its rotational and gravitational energies: 

Srot^ -§^^7 X l0-^PfRlMf„ ( 31 ) 

xigrav 

where R = R(,x 10** cm, B x 1 s, and M = x 1.4Mq. This 
deformation is pointed along the rotational axis. 

However, not the entirety of this ellipticity participates in pre¬ 
cession. The neutron star supports hydrostatic stresses only after 
the crystallisation of its crust. The part of the ellipticity that partici¬ 
pates in precession is much smaller than the value given by eq. (31) 
and can be written as (see, e.g., Munk & MacDonald 1975): 

Ccr = = 2 X lO-"p 3 oP-^RlMf„ (32) 

where p = \9p‘°°jlpgR, where g is the surface gravity, p is the 
density and p“ is shear modulus of the crust, whose fiducial value 
is lO*** dyn cm“2 

The impact of the magnetic field can be calculated in a similar 

way: 

Smag 10-^^B\XM-^1, (33) 

•^grav 

where B = B 12 x 10*2 G. This deformation is pointed along the 
magnetic moment. We discuss the possible origin of stellar non¬ 
sphericity in Section 6. 


5.2 Pulsar B1828-11 

Pulsar B1828-11 was the first to show highly periodic long-term 
variations in timing residuals that are correlated with variations in 
the pulse profile. The spectrum of timing residuals consists of three 
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Table 2. Pulsar B1828-11: Best fit parameters obtained by fitting the pulsar 
timing residuals (see Fig. 4). 


Case 

8, deg 

L. deg 

£^13 

large y- 

5 

89 

9.4 X lO-** 

small X 

89 

5 

5.4 X 10-’ 


harmonics with periods 1000, 500 and 250 days. If the period of 
free precession Pjs (see eq. 23) equals the largest period in spec¬ 
trum, ellipticity equals e =» 4.69 x lO^**. 

We will use the following timing parameters for pulsar B1828- 
11; P® = 0.405 s, = 60. The values with the superscript ‘0’ 
denote the best fit parameters, and the values without this super¬ 
script denote the observations. The characteristic timescales of the 
pulsar evolution are: 

r = 6.75 X lO'^if s (34) 

Tanom = 3.49 X lO’^fc, ' S, (35) 

where t is spindown timescale given by eq. (16) and ^ = 1 -l- sin^ a 
is a number between 1 and 2 that reflects the dependence of MHD 
pulsar spindown on the inclination angle (see eq. 5). 

This pulsar also has a well-measured second period derivative 
= -1.70 X 10“^^ and the changes in the pulse profile are 
correlated with the changes in P and P. As seen in Fig. 4, the period 
timing residual, AP = P - P^ - P^t — P^f 12, is periodic and has 
an amplitude of 1.65 ns, and the residual for the period derivative, 
AP = P - P® - P®r, is also periodic and has a magnitude of about 
0.19 X 10“*^. In order to compare the outcome of our models to 
the observations, i.e., compare the simulated values of AP and AP 
to the observations in Fig. 4, we will use the same procedure as 
was used to determine the observed values of residuals (Lyne et al. 
2013): every 50 days (indicated by dots in the figure), we average 
the residual over the time interval of 100 days. 

Stairs et al. (2000) concluded that the most probable expla¬ 
nation for these variations is the precession of the neutron star. 
Link & Epstein (2001) presented a calculation, where they at¬ 
tempted to find the best-fitting pulsar geometry under the vacuum 
approximation, that matches the residuals, AP, and the changes in 
the emission profile. We now look for a solution under a more real¬ 
istic assumption of the plasma-filled magnetosphere. 

To reproduce the observational data for pulsar B1828-11, we 
set Tprcc = ±0.218 years. The timescale can formally be negative 
because of the negativity of s; for the comparison of time scales, 
what matters is its absolute value, which corresponds to the period 
of variation of about 500 days - the leading harmonic in the spec¬ 
trum - and e = ±9.44 x 10“®. One can also show that for a biaxial 
star, the timing residuals do not depend on the sign of the ellipticity. 

We numerically solve equations (20)-(22) for different initial 
values of a, 6 and;y. After that, we find the mean value of the period 
and its first two derivatives, and compute the residuals: 

SP = P-{P}- {P)t - (P)d 12, (36) 

6P = P-{P)-{P)t. (37) 

Angles 6 , x and o are set to their initial values (which are 
denoted with the subscript ‘0’) at time t = 0 that corresponds to 
50,300 MJD. Euler’s equations are solved in t from -3 years to 
3 years to mimic the observed data range. The mean values are 
calculated over the entire time range. 

Our best fit gives x = ^9 deg and 9 = 6 deg and is presented 
in Eig. 4 and in Table 2. The solution implies that the evolution of 
a neutron star is close to free precession with small perturbations 


caused by the magnetospheric torque. Inclination angle a varies 
from;(f - 9 = S3 deg tox + 6 = 95 deg and its initial value simply 
sets the initial phase of the evolution. As we will discuss below, 
a symmetry exists between the values of 9 and x- This symmetry 
implies the existence of a “mirror” solution with the values of 9 
and X swapped (see Table 2). To distinguish between these two 
solutions, we refer to them as the “larger” and “smallsolutions. 
For brevity, we do not show the small qf solution in Figure 4. 

In order to reproduce the magnitude of the observed residuals, 
in our numerical solution we need to choose the angle 9 to be small 
and the angles large, or vice versa: as we discuss below, this is re¬ 
quired for the residuals to have such a small amplitude as observed, 
as well as the two-peak structure of the P residual. 


5.3 Analytical solution for B1828-11 

Our numerical solution for pulsar B1828-11 implies that the mo¬ 
tion of a slightly non-spherical neutron star is essentially a free 
precession with a small perturbation caused by the magnetospheric 
torques.* 


5.3.1 Free precession 

For freely precessing body, Euler’s equations of motion take the 
following form: 

L ± a X L = 0. (38) 


These equations can be solved analytically in terms of Jaco¬ 
bian elliptic functions cn, sn and dn (Landau & Lifshitz 1976): 


Li 1 10.(1 = sin So cn(a)pt, ktan^ 9(,), 

L 2 IIQ .0 = sinSo(l ± L)'^^ sn((jjpf, Ltan^ 60 ), 
L 2 IIO 0 = cos So dn(tUpf, fctan® So), 


where 


and 


k = 


Sl2 


h(h — h) _ gi2(l + gp) _ 

hih-h) gi3 - gl2 gi3 - gl2 


fi^LcOsSo £i 3 floCOsS 


h{\+kfF 


(39) 

(40) 

(41) 

(42) 

(43) 

(44) 


(1 ±fc)‘/2 ■ 

The precession period is equal 

_ P 2E(;r/2,Ltan^So) 

gi3 COS So n 

where F{^, m) is the Legendre elliptic integral of the first kind. 

In the case of £12 = 0, that is of an axisymmetric star, the 
solution takes a much simpler form: 

Li //flo = sin So cos(Wpf), (45) 

Lo/ffio = sinSo sin(Wpf), (46) 

Lo/fflo = cos So, (47) 

where Wp = gnOo cos So and Oq is the initial rotational frequency. 
In this case, the precession period is simply 

P 


(biaxial star) 


(48) 


£13 COS So 

Now we have an exact solution f2o(0 for free precession and can 


* So long as we seek a solution for the time interval much smaller than the 
alignment time. 
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t, years t, years 

Figure 4. Numerical solution of Euler’s equations (blue squares) giving our best-fit for PSR B1828-11, = 89 deg, 6 = 6 deg, and the observed residuals 
(red dots). [Left panel:] Period residual. The numerical solution is in good agreement with observations, reproducing both the characteristic period and the 
amplitude. [Right panel:] Period derivative residual. Our numerical solution reproduces the two-bump structure well. 


use perturbation theory to find the solution for a nonzero magneto- 
spheric torque. 


5.3.2 Full perturbative solution for a biaxial star 
We find the exact solution for the equation 

L + nxL = K [no(f)]. (49) 


in Appendix A. Here we summarise the most important results. 

For an axisymmetric (biaxial) MHD pulsar, the residual in P 
takes the following form (see Appendix A4 for a derivation): 


SP = (00,x)lco^ + g(do,x) cos 2wpt] 

where is the observed value of P and 

^ sin 29 sin 2 x 

4 - 2 cos^ 0 cos^ X - sin^ 0 sin^ ^ 
g{e,X) = tan0tan;i^/4. 


(50) 

(51) 

(52) 


The amplitude of the residual is proportional to sin 20 sin 2x, 
and the residual does not depend on the anomalous torque. If we 
have residuals that are much smaller than (which is usually the 
case), the best-fit values of 0 and;^^ are pushed close either to 0 or to 
90 deg. As the observed residuals for B1828-11 are much smaller 
than and we observe the two-hump structure (and thus have 
g(d,X) ~ l)i one of the angles is close to 0, and the other one close 
to 90 deg. 

One can see that expressions (51) and (52) are symmetric in 
respect to the substitution 9 x- This means that we cannot dis¬ 
tinguish between the cases of small 0 and small x, which, as we 
mentioned previously, we refer to as the cases of large x and small 
X, respectively. For small 0, the stellar deformation points in the di¬ 
rection of stellar rotation and it is natural to assume that it is caused 
by rotation. If so, it should be on the order of the analytic expec¬ 
tation, (eq. 32). On the other hand, for small it is natural to 
assume that the deformation is caused by the magnetic field. If so, 
we would expect the ellipticity to be on the order of eamg (eq- 33). 

The special case of a biaxial star corresponds to k = 0 (see 
eq. 42). For k = 0 and a small 0 value, we have cos0 ~ 1 and 
Tprec = Pis (see eq. 48). This implies s = 9.4 x 10“'*. For k = 0 and 
small we find cos 0 ~ (;r/2 - 0) and Tp^c = Pls{7tl2 - 0). For 0 = 
89 deg, this implies e = 5.4 x 10“^. These results are summarised 



Figure 5. Method of obtaining geometrical parameters of processing neu¬ 
tron stars from the observations. One needs to measure the residual values 
at three extremum points (Aa, Ab, Ac) and then calculate fix, 6) and g(x, 6) 
using equations (53)-(55). After that, x and 6 are implicitly set by the values 
of / and g. 


in Tab. 2. We discuss the possible origin of the ellipticity in Sec. 6. 
Note that for our numerical solution shown in Fig. 4 we chose a 
small 0 and large For brevity, we do not show the opposite case. 

Using the analytic model we have developed, we now can es¬ 
timate directly from the observations the values of angles 0 and x, 
or, equivalently, the values of the functions f(9,x) and g{9,x) (see 
eqs. 51 and 52). For this, we only needs to measure the following 
three extremum points, as shown in Fig. 5: the global maximum Aa 
and minimum Ab as well as the local minimum Ac of the residual. 


Aa 

Ab 

Ac 


f(So,x) 


8i8o,x) + 


1 


8g(0o,T) 

-m,x)n + g(Oo,x)i 

f(do,x)[i - gi8o,x)l 


(53) 

(54) 

(55) 


invert these to obtain the values of / and g. Note that this system 
of equations is over-constrained, and this offers a self-consistency 
check. 


5.3.3 Full perturbative solution for a triaxial star 

If the star is not axisymmetric, the solution has a much more com¬ 
plex form (but still does not depend of the anomalous torque, see 
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t, years from now 


t, years from now 


t, years from now 


Figure 6. Blue line shows our best-fit numerical solution of Euler’s equations for the Crab pulsar in the large x is shown (see Fig. 7 for the small case). 
This red line indicates the time range over which timing residuals have been observed and the black dot indicates the present time. Black dashed lines show the 
observationally-inferred values. The pulsar parameters at f = 0 are: rro = 70°, Qq = 2°,xo = 68.05°. [Left panel:] The rate of change of the inclination angle, 
a. The solution shows an order of magnitude agreement with observations. [Central panel:] The braking index /tbr vs time. [Right panel:] The inclination 
angle a vs time. Over the observed time range, a increases almost linearly with time, in agreement with observations. 


Appendix A5 for a derivation): 

6 P ^ -P°'^’‘f(eo,x)n0o,x,t) 


(56) 


with 


P = cn{ci)pt)dn{ci)pt) + cn{2(jjpt) 
+ (dn(2wpt) - 1) 


2g(eo,^) ^ 

(1 + dn(2aipt)) 

gido,x) , gido^xL' 


(57) 


(1 + dn(2(j>„r)) 8(1 + cn(2a)pl)) 


Here we for simplicity do not write the second argument of Jaco¬ 
bian elliptic functions which always equals ktan^ Og. The second 
ellipticity does not qualitatively change the behaviour of solution 
and only influences the precession period (see eq. 44). 


5.4 Crab pulsar 

The emission from the Crab pulsar (pulsar B0531-I-21) was ob¬ 
served since 1969 in radio, optical, x-ray and gamma-ray wave¬ 
bands. The radio pulse profile of this pulsar consists of two com¬ 
ponents separated by 145° in phase: the main pulse and the inter¬ 
pulse. During high-precision daily observations performed since 
year 1991 at Jodrell Bank Observatory, Lyne et al. (2013) found 
a steady increase in separation between the main pulse and the in¬ 
terpulse, at 0.62° ± 0.03° per century. 

Because radio and gamma-ray pulses of the Crab pulsar ar¬ 
rive at the same pulse phase, it is thought that the radio emission 
from the Crab pulsar comes from the same magnetospheric region 
as the gamma-ray emission. Using this fact, as well as gamma-ray 
emission models, Lyne et al. (2013) concluded that the increase in 
the separation between the main pulse and the interpulse indicates 
an increase of the inclination angle at approximately the same rate. 
The increase of the inclination angle also naturally leads to a brak¬ 
ing index value less than 3, in agreement with observations. How¬ 
ever, for a spherically symmetric star, MHD models predict a de¬ 
crease of the inclination angle for a spherical star (see Sect. 4), and 
a braking index that is always larger than 3. Philippov et al. (2014) 
suggested that the observed increase of the inclination angle could 
be due to the precession of the neutron star caused by stellar non¬ 
sphericity. 

Later, Lyne et al. (2015) revealed the observational data for 
the Crab pulsar over a 45-year time range. Although the timing 
residuals show some oscillatory behaviour, they are corrupted by 


a large number of glitches. This makes it impossible to get an ac¬ 
curate precession model fit. Nevertheless, between the subsequent 
glitches Crab pulsar behaves as if it were precessing. 

The observed increase in the inclination angle of the Crab 
pulsar can be explained with a precession model. Just as for PSR 
B1828-11, the timing residuals residuals in P and P (Lyne et al. 
2015) imply that either S x; 0 and x ~ 90 deg or s; 0 and 
0 x 90 deg. 

Figure 6 shows our best fit obtained for large x case and 
a = 70°. One can see that for the past 45 years the inclination of 
the Crab pulsar has increased almost linearly. On larger timescale, 
the inclination angle is oscillating and its mean value is almost con¬ 
stant. Our fit produces values of v, v and v close to the observed 
values. We summarise the inferred pulsar geometry in Table 3. The 
precession period in this solution is roughly Tprec = 126 years, but 
there is a strong degeneracy between values of 6 and Tprec- Making 
precession period Tpiec larger (which is allowed by the data) and the 
wobble angle smaller, one could obtain sufficiently good fits. 

Fig. 7 shows the solution for precession model in the case of 
small The inferred geometrical parameters are listed in Table 3. 
The solution reproduces the observed data for the past 25 years, 
but unlike the case of larger, the inclination angle oscillates on top 
of a global alignment trend. There is the same degeneracy between 
0 and Tprec. So, ~ 100 years is roughly the smallest period which 
leads to the good agreement with observations. 

For isolated pulsars it is common to consider a dimensionless 
quantity called the braking index: 

on 

n.-—. (58) 

A spherical pulsar rotating in vacuum has 

„VAC = 3 + 2tan-V>3, (59) 


and a spherical MHD pulsar has 


n 


MHD 

br 


= 3-1-2 


sin^ a cod' a 
(1-1- sW ad 


3 < n 


MHD 

br 


< 3.25. 


(60) 


Although njjj ^ 3 for both models, observations show values not 
only less than 3 (e.g., the Crab pulsar has nbr ~ 2.5) but even neg¬ 
ative values. In fact, braking indices span the range from -10® to 
-1-10®. If the pulsar experiences precession, it can have an arbitrary 
value of braking index both larger and smaller than 3, and can also 
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small X small x small x 





Figure 7. Similai' to Fig. 6, but for the small x case. The parameters of this solution at r = 0 are: a = 60°, 6 = 59.1°, x = 1°- The solution shows order of 
magnitude agreement with observations. Note that for^ <s: 1, i.e., when the magnetic and stellar symmetry axes are close to each other, a shows a secularly 
decreasing trend, as opposed to the large case, shown in Fig. 6, where a oscillates about a constant value. 



Figure 8. Comparison of our model results and the observed frequency residual for Crab pulsar. [Left panel:] The frequency residuals vs time in our numerical 
simulations are shown with the red line for the case of large x Ihe blue line for the case of small x- [Right panel:] The observed frequency residual vs 
time. Despite a large number of glitches in the data, our model shows an order of magnitude agreement with the observations. Note that the variations in 
the frequency residual occur not on the precession period but on a much larger time scale set by the fourth-order frequency derivative. In order to observe 
variations caused by precession (as in the case of PSR B1828-11), one needs to observe Crab pulsar for at least several precession periods. 
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Figure 9. The compaiison of our Crab modelling results to the observations of the frequency derivative residual. [Left panel] The red line shows the simulated 
residuals for the largecase and the blue line for the small^ case. [Right panel] Observed timing residuals. As in Fig. 8, our modelling results show an order 
of magnitude agreement with the observed values, as expected given the large number of Crab’s glitches that are not included into our modelling. 
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Table 3. Crab pulsar: Best fit parameters for large and small ^ cases. The second and third frequency derivatives agree with observations to within an order 
of magnitude, as expected given the large number of Crab’s glitches. The best-fit values of the inclination angle a are agreement with independent estimates 
(Lyne et al. 2013). Our solution implies a precession period of about 100 years. Deviations from the near-linear increase in time of the inclination angle could 
be observed on shorter timescales, perhaps as short as ~ 20 years. 


Case 

6, deg 

Y.deg 

a, deg 

Tprec, years 

£ 

V, 10^^'' Hz/s^ 

V, 10-'™ Hz/s ' 

observations 

- 

- 

- 

- 

- 

1.12 

-2.73 

larger 

2 

68.05 

70 

126.7 

8.3 X 10-‘2 

1.15 

-0.8 

small X 

59.1 

1 

60 

108.5 

1.11 X 10-" 

1.27 

0.8 


be negative. As we show in Figures 6 and 7, the observed braking 
index of the Crab pulsar can be easily reproduced by our model. 

We note that our precessing solutions show that the near-linear 
increase of the inclination angle of Crab pulsar will be inevitably 
followed by its decrease, and the minimal timescale for this switch 
is ~ 100 years. It is possible that deviations from the linear increase 
could be detectable on much shorter timescales, perhaps as short as 
~ 20 years. 

For the Crab pulsar, it is quite difficult to fit the timing residu¬ 
als data due to the presence of glitches that are not included into our 
modelling. Nevertheless, we can compare the characteristic ampli¬ 
tudes of simulated to observed timing residuals, as shown in Figs. 8 
and 9. In both small and large x cases, our modelling results and 
observations agree to within an order of magnitude. 

We conclude that precession is a plausible explanation of ob¬ 
served behaviour of Crab pulsar, quantitatively reproducing the 
range of observed frequency derivatives, the braking index, incli¬ 
nation angle derivative, and the timing residuals. 


6 DISCUSSION 

In this Section, we discuss and compare our results for the Crab pul¬ 
sar and pulsar B1828-11, and discuss astrophysical implications. 


6.1 Vacuum vs. Plasma-filled Magnetospheres 


We start by discussing the differences between the vacuum and 
MHD models in terms of their effect on the timing residuals. As 
we showed in Sec. 5.3, the anomalous torque has no effect on the 
timing residuals. This means that the uncertainties in its estimate 
do not affect our results. From the point of view of timing residu¬ 
als, the difference between vacuum and MHD models is evident in 
the denominator of equation (51): 


MHD ^ sin 200 sin 2^ 

^ A-A(e^,X)' 


( 61 ) 


which for the vacuum case takes the following form: 


VAC ^ sin 26 )osin 2 ^ 

2-A(0o,^)’ 

where A(6,x) = 2cos^ Bcos^x + sin^ 0sin^;^'. Since the observed 
timing residuals, given by eq. (50), are observationally inferred to 
be very small, i.e., 6P/P « 1, the / factors given by eqs. (61) 
and (62) are also « 1. To get a good fit, the values of 9 and x 
must be either near zero 0 or 90 degrees. Importantly, the two-peak 
structure, such as seen in the residuals of the pulsar B1828-11, is 
reproduced only for g = tan;ytan6o/4 ~ 1 (see eq. 52). So, either 
9 is close to 0 and^ is close to 90 degrees or vice versa. Thus, the 
angular dependence in the denominators of equations (61) and (62) 
is negligible, i.e., A <K 1, and 


yMHD _ Q 5 


(63) 


The twice as small residuals implied by plasma-filled mod¬ 
els lead to a larger space of allowed solutions (Sec. 5.3) and allow 
values of X S to be further away from the extreme values (e.g., 
0 or 90 degrees) than in the vacuum models. For instance, previ¬ 
ous studies of pulsar B1828-11 in the context of the vacuum model 
(Jones & Andersson 2001; Link & Epstein 2001) found 9 < \ de¬ 
gree while our plasma-filled model gives 9 = 5 degrees. 

The effect of stellar non-sphericity on the long-term evolution 
of pulsar parameters depends on the magnitude of the x angle be¬ 
tween the magnetic moment and the principle axis of the star. In the 
case of large ;y, the angles 9 and a oscillate at the precession period 
and their averages show no long-term trends, as seen in Fig. 6. 

In the case of small x, the values of a and x also oscillate 
on the precession timescale but undergo a secular decrease over 
longer timescales, as seen in the right panel of Fig. 3: both vac¬ 
uum and MHD pulsars evolve toward alignment. Thus, if most pul¬ 
sars are in the small x case, non-sphericity effects will have very 
few implications for the pulsar statistics, which favors an align¬ 
ment trends on the timescale of 10*-10^ years in agreement with 
spherical pulsar models (Tauris & Manchester 1998; Young et al. 
2010; Gullon et al. 2014). The vacuum pulsars approach the align¬ 
ment exponentially fast and the MHD pulsars do so much more 
slowly, as a power-law in time. These long-term evolutionary trends 
are quantitatively similar to those for spherically-symmetric pul¬ 
sars (Philippov et al. 2014). However, unlike the spherical pulsars, 
which evolve toward complete alignment (a = 0), non-sphericity 
causes the inclination angle to stabilise at a finite value, ctstab > 0. 
The precise value of ttstab can differ between vacuum and MHD 
models: for instance, for the parameters chosen in Fig. 3, an MHD 
pulsar stabilises at ttjf™ = 3.46 deg and a vacuum pulsar at a much 
smaller value, = 0.35 deg. Note that since vacuum pulsars 
align exponentially fast, they can reach stabilization after 10* years, 
which is much shorter than the pulsar lifetime. In contrast, stabili¬ 
sation can take ~ 10* years for MHD pulsars, and it is unclear how 
many MHD pulsars reach stabilisation in their lifetime. 


6.2 Crab pulsar vs. BI828-1I 

We demonstrated that the increase in time of Crab pulsar’s inclina¬ 
tion angle could be easily explained if this pulsar experiences pre¬ 
cession, confirming our earlier expectations (Philippov et al. 2014). 
A precession period, Tp^c ^ 100 years, and rotational period, 
P = 0.033 s, imply a rather small stellar ellipticity scrab ^ 10“" 
(see Tab. 3). Indeed, this value is much smaller than the minimum 
allowed value of ellipticity for the pulsar B1828-11, ei 828 ~ 10“*. 

When this work was completed, a preprint by Zanazzi & Lai 
(2015) was posted on the archives. They considered a precession 
model for the Crab pulsar in the framework of vacuum magneto¬ 
sphere and obtained e = 4 x 10“", slightly higher than our value. 
Their geometrical parameters are more extreme, as one would ex¬ 
pect from a vacuum model: x = 0.15 deg (vs <K 1 deg in MHD 
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models) in the small x case and 6 = 0.1 deg (« 2 deg) in the large 
^ case. 

For the Crab pulsar, the MHD spindown formula (11) implies 
a surface magnetic field S 12 = 2.33 for a characteristic pulsar obliq¬ 
uity angle ff = 60 deg (the magnetodipole formula gives a slightly 
larger value B 12 = 3.78). The estimate (33) for magnetically- 
induced ellipticity gives s = 0.6 x 10“**, i.e., very close to the 
inferred value, so Crab’s ellipticity could be accounted for by a 
magnetically-induced deformation of its crust. We note that the we 
can only place a rough lower limit on the precession period of the 
Crab pulsar, and it could be significantly longer than 100 years. 
Long-term observations in the next several decades might help to 
break this degeneracy and measure the precession period if it is not 
much longer than a century. 

The pulsar B1828-11 has approximately the same surface 
magnetic field, but its inferred ellipticity is 3 orders of magnitude 
larger, e — 10“*. One way to reconcile this difference from the Crab 
pulsar is to postulate strong (~ 10*"* G) internal magnetic fields in 
the pulsar B1828-11. Another possible explanation for the differ¬ 
ence in the values of ellipticities can be found in the framework 
of rotational ellipticity model. Equation (31) gives very large ellip¬ 
ticity values. However, as Jones & Andersson (2001) pointed out, 
only a small fraction of this deformation participates in the pre¬ 
cession due to the pinning of magnetic vortexes to the crustal lat¬ 
tice. Crab is younger, has a lot of glitches, so, the pinning fraction 
is large and the effective ellipticity small. Since PSR B1828-11 is 
older, it is conceivable that it has a smaller pinning fraction and a 
larger effective ellipticity. In this context, it is interesting that all 
known pulsars that experience timing noise on the timescale of 
several years, have ages larger than that of B1828-11 (Lyne et al. 
2010). Thus, we may expect that the older the pulsar, the larger the 
effective ellipticity that participates in precession. Dissipative pro¬ 
cesses in the stellar interior such as crust-core coupling may also 
affect the long-term evolution of the neutron star (Barsukov et al. 
2014). 

Pulsar B1828-11 does not show an observed interpulse, and 
this might seem at odds with our precessing solution, in which a is 
very close to 90 degrees. However, due to the uncertainties in the 
radio emission mechanism and the plasma distribution near the po¬ 
lar cap of nearly orthogonal rotators, it is not clear that an interpulse 
will necessarily be visible. It is possible that magnetospheric syn¬ 
chrotron absorption (Beskin & Philippov 2012) can significantly 
damp the interpulse emission. More detailed investigation is needed 
to produce reliable models of the observed radio profiles of this 
pulsar and independently confirm the geometry inferred from pre¬ 
cession models. 


7 CONCLUSIONS 

In this work, we presented the analysis of pulsar evolution that self- 
consistently takes into account magnetospheric plasma effects, as 
measured in 3D MHD simulations of pulsar magnetospheres. We 
showed that the MHD effects expand the allowed parameter space 
of pulsar precession solutions and thereby allow for less extreme 
solutions than the vacuum models. To facilitate the interpretation 
of pulsar timing residual curves, we developed an analytic model 
that converts the observed residuals into the pulsar geometrical pa¬ 
rameters. 

We applied this model to the timing residuals of PSR B1828- 
11 and obtained its parameters: e.g., stellar ellipticity and the ori¬ 
entation of the stellar ellipsoid, magnetic and rotational axes. With 


the best-fit parameters, both the period and period derivative timing 
residual curves are well-reproduced by our model. Our best-fitting 
values for the pulsar parameters are less extreme than in the vac¬ 
uum models. 

In the context of spherically-symmetric stars, it is expected 
that the magnetic and rotational axes of the pulsars evolve toward 
each other as the pulsar ages. Recently a puzzling observation of an 
opposite trend was reported for the Crab pulsar (Lyne et al. 2013). 
In this work, we showed that this surprising trend can be explained 
by the precession of a neutron star, in agreement with our earlier 
expectations (Philippov et al. 2014; see Zanazzi & Lai 2015 for a 
consideration of this effect in the vacuum approximation). The data 
is consistent with a precession period of the Crab pulsar that is as 
short as < 100 years. If the true precession period is not exceed¬ 
ingly larger than that, it could be measured in the next decades, 
warranting a continuous monitoring of the Crab pulsar. 
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APPENDIX A: ANALYTICAL SOLUTION 

In this section we describe analytical solution for pulsar precession. 
The exact solution could be treated as small perturbations of free 
precession motion. We describe free precession in Section Al, per¬ 
turbation caused by anomalous torque are discussed in Section A2, 
effects of regular alignment torque are considered in Section A3, 
two-peak structure in residuals are obtained in Section A4 and we 
finally discuss the effects of the second ellipticity in Section A5. 


Al Free precession 

For freely precessing body Euler’s equations of motion take the 
following form: 

L + ilxL = 0. (Al) 


These equations can be solved analytically 
(Landau & Lifshitz (1976)) in terms of Jacobian elliptic functions: 


LiIIQq 


L 3 IIQ .0 


sin 60 X cn(£jpf, k tan^ Oq), 
sin 9o dl + k x snftUpt, k tan^ 9o), 
cos 00 X dnftUpf, k tan^ 0o), 


(A2) 

(A3) 

(A4) 


where A: - h{h—h)lhi.h—h) - £i 2 (l+£i 3 )/(£i 3 “£i 2 ) ~ £i 2 /(£i 3 “ 
£ 12 ) and (Up = £i3Lcos0o/ 73 Vl + ^ ~ fiisOo cos 0/ Vl +k. The 
precession period is equal 


^prec — 


P 2F(nl2,ktardeo) 
£13 cos 00 n 


(A5) 


where F(<p, m) is the Legendre elliptic integral of the first kind. 

In the case of £12 = 0, that is of an axisymmetric star, the 
solution takes on a much simpler form: 


Li 110.(2 

= sin 00 X cos{(Opt), 

(A 6 ) 

L 2 lIO (2 

= sin 00 X sin(o)pf). 

(A7) 

L 2 IIO 0 

= cos 00, 

(A 8 ) 


where Wp = fioOo cos 0 o and flo is the initial rotational frequency. 
In this case, the precession period is simply 

fprec = P/si 3 COS 9q. (axisymmetric star) (A9) 


A2 Effects of the anomalous torque 


In reality the precession of radio pulsar is not free. In order to cor¬ 
rectly describe it one have to solve the following equation: 


L + SlxL = K = 

“ Aaiigned 


(AlO) 


11 n k-i sm a cos a Fix u 

-2 - v cos a— H- 

12 /j SlR/c |I2Xyu| 


where we already use MHD torques with ko = ki = k 2 = \. 

After we normalize this equation by IFlo and set L = I(Sl + 
£1212262 -t- £ 1212363 ), we obtain the following expression: 


1 d >3 £12 CO2 612 

0} + - "I" X ^3 + - COt,OJ X ^2 ~ 

Tprec I 2 o £13 I 2 o ' £13 

(All) 


(o n 

- -2-h cos O' — 

TO) 


ad (OX ft 

-I-sm O'cos a---, 

^anom 1^ ^ /^l 


where m,- = n,/I2o, Vec = t = ^i^oZ-S'^gned '^anom = 

(I2„2?/c)/n„/*3X«jg„,,. 

In this section we consider £12 = 0 and Tp^c « « t. 


This approximation is well-satisfied for pulsars and we can treat 
torque effects as perturbations of free precession motion. 

Using these assumptions and the relation costt = Wi sin;^f - 1 - 
(02 cos;^' we can write the equations for perturbations (5a), in the 
following way: 

0 ) 3 ( 50)1 -I- 0 ) 100)3 
fprec 

ad 

-( 0)1 sm;^^ -I- (02 cos;^')o )2 cos;^', (A12) 

fanom 

(02,6(02 + (026(02 

f piec 

id 

-(tiJi sin;^" -I- 0)3 cos^)(o)3 sin;^^ - o)i cos^), (A 13 ) 

^anom 

ad 

-(tiJi sin;^'-I-0)3 cos^)o)2 sin;^f. (A 14 ) 

^anom 

Here we neglect the effects of t and the first two terms in brackets 
in LHS of equation (Al 1) as they are much smaller than (o. 

Now we set o), to be equal to their unperturbed values and 
solve the system of linear differential equations (A12)-(A14). This 
is the system of linear differential equations with constant coeffi¬ 
cients, so it is straightforward to find a solution: 


00)1 + 


00)2 - 


00)3 = 


0 o)i = — (/(Q-t t/'iO)p? sin(Upf — ^/(2 cos o)pf— 

- i /(2 sin^ (Op cos (Op, 

00)2 = - Ipi(op cos (op + p 2 sin (op cos^ (op 
00)3 tan 00 cos (op -l- p 2 tan 9q cos^ ojp, 
where (o„ = cos 0 o/Tprec and 


do = 
Pi = 
p2 = 
lp2 


■ sm;^fcos;^'cos0o, 

2 2 2 sin 00 - sin^ 0o 

cos X sin 00 - sin x 


4 cos^ 00 


‘prec .2 ■ r 

-sm ;^^sm 0 o/ 2 . 


Tprec sin^ 0o 


(A 15) 
(A16) 
(A17) 

(A 18) 
(A19) 
(A20) 
(A21) 


Tanom 4 COS^ 

It is important to notice that the amplitude of (o remains unper¬ 
turbed: 0 ) 100 ) 1 -I- 0 ) 200 ) 2 - 1 - 0 ) 300)3 = 0. This is because the anomalous 
torque is perpendicular to 12 and thus cannot change its amplitude. 

On the other hand, anomalous torque causes the changes in 0 
since 00)3 = - sin 9q69. 

cos (0„t 


69 = ■ 


-iPo + p 2 cos (op). 


(A22) 


The maximum deviation is given by 


max |0 - 0o| 


Tprec sin 2x + sin^^ tan 0 o 


(A23) 


A3 Effects of the alignment torque 

If one considers in equation (All) the contribution of alignment 
torque, it is possible to complete exactly the same procedure as in 
Section A2. The solution in this case is quite different: 

00)1 = (p2sm2(Opt - (/) 2 (Optcos(Opt + (p 4 (Opd sin (Opt, (A24) 
00)2 = (j>0 + (f>l cos (Opt — (p2 cos 2o)pf — (^ 30 )pf sin (Opt— 

- (p4(0^pd cos (Opt, (A25) 

00)3 = {(p 2 sin (Opt - (putOpf) cot 0o, (A26) 
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- numerical 

solution 

- anomalous + 

alignment 

-anomalous 

alignment 


Figure Al. The compaiison of exact numerical solution with analytical asymptotic for^ = 88.5 deg, Oq = 6 deg. The figure shows an excellent agreement of 
numerical solution and observations. This allows us to use simple analytical expressions (56)-(57) instead of numerical procedures. 


where 


<^o 

<^2 

03 


2'prec 2 — tcin^ 

---sin;(-cos^, 

T 2 

Tprec sW^tan^o 
r 2 

^prec . 2 ^ 

-sin;^ cos;^ tan 

T 

^prec 7 + COS 2y 

- -tan 00 , 

T 4 


(A27) 

(A28) 

(A29) 

(A30) 


04 - 


fprec 3 — COS 2^ 


tan 00 - 


(A31) 


If we combine the effects of both alignment and anomalous 
torque we obtain the total perturbation. Fig. Al compares exact nu¬ 
merical solution with total perturbation as well as the perturbations 
caused by alignment and anomalous torques separately. 

The alignment torque causes o) to change. The perturbation to 
the frequency takes the following form: 


6(l) — -t- ix)q6(02 + — 


(A32) 


= —^Wpr(3 - cos ly + 1 tan^ 0o - cos 2y tan^ 0o) cos 0o/2-l- 

T 

H—^ sin 2x sin 0o(sin (jjp +1/8 tan 9q tsax sin laip). 


Here we can see spin down as well as the two harmonics that are 
actually presented in observations. 

Summing up, the observed residual curves are the first-order 
perturbation over free precession motion caused by the alignment 
torque. The first-order residual has up to two harmonics. 


A4 Two-peak structure 


The frequency residual 6a) could be found in much more simple 
way. In order to find it we multiply equation (Al 1) by <u: 


£U 3 d >3 + £{ 2 /^ 130 ) 26)2 n3^(l + sin^ cr) 




(A33) 


Puting zero-order solution for £12 = 0 in RHS of this equation 
gives the following equation: 


6) = —(2 - cos^ 00 cos^;^' - 1/2 sin^ 9q sin^;^')+ (A34) 

T 

1 sin 200 sin 2x 

+ - - -(cosWpf + cos2wpttan0o tan;^^/4). (A35) 

The first term of this expression is the observed spindown. 
Having that we can rewrite this equation in the following way: 


/ sin 200 sin 2;^f(cos ajp + cos 2o)pt tan 0o tan;^^/4) \ 

P°’’* \ 4-2 cos^ 00 cos^x ~ sin^ 0o sin^;^^ / 

(A36) 


Finally, the residual in P now takes a form 


6P = -P°^^f(9o,x)(cos 0 )p + g(9o,x) cos 2a)p) 

(A37) 

with 


sin 200 sin 2^ 

fiSo,x) ~ ^ , . 2 „. 2 ’ 

4-2 cos^ ^0 cos^;^ - sin 6o sin x 

(A38) 

g(.6,x) = tan0otany-/4. 

(A39) 

The amplitude of residual is proportional to sin 20o sin 2x- If 
we have residuals which are much smaller than (which is usu¬ 

ally the case), 0o and x should be close either to 0 or to 90 deg. 
As observed residuals are much smaller than and we observe 


two-bump structure (so, g(9Q,x) ~ 1), one of the angles should be 
very close to 0, and another one - to 90 deg. 


AS Effects of the second ellipticity 

For £ 12 , £23 « 1, one can show that 

£ 12 . ~ ^ 

£13 k + 1 

In this case, zero-order solution implies 


0)36)3 + £{2/£230)26)2 = 0 . 


(A40) 


(A41) 


The solution for perturbations takes now very similar form: 

6P = -P°'^y(9o,x)P(()o,X, t) (A42) 


with 


P = cn{o)pt)An(a>pt) + cn(2wpf) 
+ {dn{2o)pt) - 1) 


2g(9o,T) , 

(1 + dn(2tUpf)) 

g(So,x) , gWo^xP' 


(A43) 


(1 + dn( 2 a)pf)) 8(1 + cn( 2 tUpf)) 


It shows that the second ellipticity does not qualitatively 
change the behaviour of solution and only influences the preces¬ 
sion period (see eq. (A5)). 


This paper has been typeset from a Tj^X/ KTgX file prepared by the 
author. 
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